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ABSTRACT 



Linear prediction has become an important tool for stationary time series analysis. 
The all-pole system model which is a product of the linear predictive approach has ap- 
plications in numerous engineering problems. This thesis develops a simple method for 
obtaining the two-dimensional all-pole system model. The lattice structures that can be 
used to implement the prediction error and synthesis filters are also shown to have an 
analogous two-dimensional counterpart. The construction of these filters is in terms of 
orthogonal Szego polynomials which can be used to solve the two-dimensional block 
Toeplitz normal equations in a recursive manner. This recursion not only leads to the 
two-dimensional (2-D) lattice structure, but also allows for expansion of the filter order 
without resolving the normal equations. Several examples are presented using the two- 
dimensional linear prediction results for spectral estimation and signal synthesis. 



TABLE OF CONTENTS 



I. INTRODUCTION 1 

A. LINEAR PREDICTION 1 

B. ONE-DIMENSIONAL LINEAR PREDICTION 2 

1. All-Pole Model 2 

2. Normal Equations and Orthogonality 5 

3. Applications of the Autoregressive Model 7 

C. NOTATIONAL CONVENTIONS 7 

1. Matrix Notation 7 

2. Causal Samples 7 

3. Polynomial Notation 8 

II. TWO-DIMENSIONAL LINEAR PREDICTION 10 

A. PROBLEM DEFINITION 10 

B. TWO-DIMENSIONAL ORTHOGONALITY 12 

1. Description 12 

2. Derivation 12 

3. Two-Dimensional Normal Equations 14 

C. SZEGO POLYNOMIALS 17 

1. Levinson Algorithm in Terms of Szego Polynomials 18 

2. 2-D Extension 20 

D. LATTICE IMPLEMENTATION 24 

1. Analysis Lattice 24 

2. Synthesis Lattice 2S 

III. COMPUTER IMPLEMENTATION AND RESULTS 34 

A. IMPLEMENTATION 34 

1. Autocorrelation Function 34 

2. Reflection Coefficients 34 

3. Filter Coefficients 35 

B. RESULTS 36 

1. System Model 36 



2. Spectral Estimation 

3. Signal Synthesis . 



36 

39 



IV. CONCLUSIONS .< 60 

LIST OF REFERENCES 61 

INITIAL DISTRIBUTION LIST 62 



v 



LIST OF FIGURES 



Figure 1. Prediction error and system model filters • 4 

Figure 2. Sample space for N1 = N2 = 4 8 

Figure 3. Causal Data Space 9 

Figure 4. Generation of the prediction error 11 

Figure 5. A visualization of the orthogonality principle for two samples 13 

Figure 6. Variables used to form the bivariate orthogonal polynomials 21 

Figure 7. First row of the analysis lattice 25 

Figure 8. The lower order polynomials for 26 

Figure 9. Coefficient relationship between orthogonal polynomials 29 

Figure 10. One cell of a 2-D analysis lattice 30 

Figure 11. First row of the synthesis lattice 32 

Figure 12. e 00 node of the synthesis lattice 33 

Figure 13. System models for Cq. 1 36 with an impulse input 37 

Figure 14. System models for Eq. 136 with a white noise sequence input 38 

Figure 15. Generating model for spectral estimates 40 

Figure 16. Surface plot of spectral estimate 41 

Figure 17. Contour plot of spectral estimate 42 

Figure 18. Surface plot of spectral estimate 43 

Figure 19. Contour plot of spectral estimate ! 44 

Figure 20. Surface plot of spectral estimate . . . 45 

Figure 21. Contour plot of spectral estimate 46 

Figure 22. Surface plot of spectral estimate 47 

F igure 23. Contour plot of spectral estimate 43 

Figure 24. Surface plot of spectral estimate 49 

Figure 25. Contour plot of spectral estimate 50 

Figure 26. Surface plot of spectral estimate 3 j 

Figure 27. Contour plot of spectral estimate 52 

Figure 28. Surface plot of spectral estimate 53 

Figure 29. Contour plot of spectral estimate 54 

Figure 30. Surface plot of spectral estimate 55 

Figure 31. Contour plot of spectral estimate 56 



VI 



Figure 32. Surface plot of spectral estimate 57 

Figure 33. Contour plot of spectral estimate SS 

I iguic 34. Spectrum of signal generated from white noise 59 



vii 



I. INTRODUCTION 



A. LINEAR PREDICTION 

The initial goal of linear prediction is to estimate the present value of a random se- 
quence based on a linear combination of all the past values. This is equivalent to de- 
termining a linear, shift invariant, causal prediction error filter which whitens the 
random sequence. There has been a great deal of study of one-dimensional linear pre- 
diction analysis and the usefulness of the resulting filters. The purpose of this thesis is 
to extend in a straightforward manner the results of one-dimensional linear prediction 
to two-dimensional signals. 

There are several efficient methods for the solution of the normal equations that 
arise in one-dimensional linear prediction problems. The most familiar of these methods 
is the Levinson algorithm. The Levinson algorithm alleviates the need to perform a 
matrix inversion in the solution of the normal equations for the filter coefficients. The 
recursive nature of the algorithm also allows the prediction error filter to be implemented 
in a lattice structure. These two features will be of primary concern in the solution of 
the two dimensional problem. 

There are two general approaches to linear prediction. The first method, which is 
usually referred to as the spectral factorization method, begins with the knowledge of the 
power spectral density of the random process and attempts to find the model parameters 
that fit the spectrum. This works well in one-dimensional applications but cannot be 
extended to two-dimensional signals. In general, a two-dimensional spectrum is not 
factorable into two polynomials. Although progress has been made in factoring two- 
dimensional spectra in terms of the complex cepstrum and the results are theoretically 
exact, these theoretical results cannot be attained in practice. We will, therefore, focus 
our attention on a second method of attaining the prediction error filter. 

The solution method dealt with in this thesis begins with the individual samples of 
the random process and develops the prediction-error filter from these samples. This 
method is known as autoregressive model fitting or all-pole modeling. The solution will 
be developed in terms of orthogonal polynomials that will not only yield a Levinson-like 
algorithm but are easily extended to two-dimensional signals. 
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B. ONE-DIMENSIONAL LINEAR PREDICTION 
1. All-Pole Model 

We are given a finite length sequence of the past values of the random sequence 
Y(n). The problem is to estimate Y(n) based on a linear combination of these samples. 

For the simplest case we can estimate Y(n) based solely on the previous sample. The 

estimate, denoted by Y(n), is 

Y{n) = aY{n- 1) (1) 

The prediction will have some error, e{ri), which will be minimized in some sense 
by the choice of a. The least-squares minimization criterion is commonly used. This 
criterion provides a simple approach to the solution of the linear-prediction problem and 
thus the modeling problem. We proceed as follows. 

First, we redefine the estimate as 

Y{n) = -aY{n- 1) . (2) 

Now the prediction error will be 

e{n) = Y{n) - Y{n) 

= Y(n) + a Y(n — 1 ) . 

The squared error is to be minimum with respect to the coefficient a. This is 
easily accomplished using derivatives as shown below. Differentiating [e 2 («)] with re- 
spect to a gives 






o jL 

2e ' x da e " 



= 0 



or 



\_Y{n) + aY{n- 1)]^-[F(«) + aY[n- 1)] = 0 
da 



which yields 



[}'(«) + aY(n- 1)][F(« — 1)] = 0 . 



(4) 

(5) 

(6) 



Finally, solving for a 



2 



a 



Y(n)Y(n- 1) 
Y{n - l)Y(n - 1) 



( 7 ) 



Since the value of a must minimize the error, e(n), for all values of n, the value 
of a in (7) could be rewritten as 



n 

1 ) 

-a = (8) 

£n/-i)n/-i) 

/= i 

Note that (8) will select a such that £c 2 (/j) is minimum. It can be seen that the 

n 

coefficient a is a function of the autocorrelation elements of the given sample space, and 
this method is usually referred to as the autocorrelation or Yule-Walker approach. 

If the order of the predictor is now increased such that 

Y(n) = - [a x Y(n- 1) + a 2 Y(n - 2) + ... ] (9) 

the prediction error becomes 

e„ = }'(«)- f(«) = l’(») + V(/7- 1) + a 2 Y(n — 2) + .... . (10) 



If (10) is considered to be the difference equation of an all-zero filter, the 
prediction-error filter can be described by taking the z transform of (10) as 

E{z) = Y(z)A(z) (11) 

where 

A(z) = 1 + ctjZ + ctjz + — . (12) 

Now if (10) is rewritten and e(n ) is taken to be the input to the system, the es- 
timated signal generator model results 

Y(z) = E(z)D(z) = -^yZT(z) (13) 



where 
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D(z) = 



(14) 



1 

1 + (l\2 + CljZ + •••• 



It can be shown that e„ is a white noise sequence with variance £[<£]. This im- 
plies that Y(z) and, consequently, Y(n) can be generated from a white noise sequence 
as shown in Figure 1. 




This assumes that the sequence A(z) is causal and causally invertible. If this is the case, 
we can generate the error sequence from Y(n) and vice-versa by an invertible filter op- 
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eration [Ref. 1]. The signal generator model, B(z), is the best estimator in a least-squares 
sense and A(z) is the prediction-error filter. These filters can be used in a number of ways 
and the simple solution for the filter coefficients is the primary reason for the popularity 
of linear prediction in signal processing. 

2. Normal Equations and Orthogonality 

Recalling that we wish to minimize the sum of squares of the error sequence 
leads directly to an important result in linear prediction known as the orthogonality re- 
lationship between the sample space and the error. This principle states that the error 
e{n) is orthogonal to all the samples of Y(n) used in obtaining the estimate Y{n). That 
is, 



£[>(«), Y{n - /)] = 0, 1 < / < N 



(17) 



where N is the order of the predictor. This result is derived for the two-dimensional case 
in Chapter 11. It is the central idea behind the solution of the linear prediction problem 
in terms of orthogonal polynomials. It also leads directly to the normal equations in 
both one- and two-dimensional formulations. 

The derivation of the normal equations is a simple extension of the 
orthogonality relationship (17). Suppose we seek a kth-order predictor and the value 
of the prediction-error variance. We will require (k + 1) linearly independent equations 
to solve for the k prediction coefficients and the prediction-error variance. These coef- 
ficients can be derived as follows. 



e n = >’(«) ~ Y('0 

= Y(n) + a, Y(n - 1) + a 2 Y(n - 2) + ■ • • • + a k Y{n - k ) 



(IS) 



Substituting (IS) into (17) gives 

£[ }'(/*) + a, >'(« - 1) + a 2 Y{n - 2) + - , }'(« - /)] = 0 1 < i <, N (19) 

If we now define the coefficient a 0 = 1 and note that for a stationary process the 
autocorrelation, £(„r), is given by £[T(A-)T(/)D = £*-/ , (19) may now be written as 



v 



) a k Y{n-k), Y(n - /) 



= 0 



/<= o 



(20) 



5 



( 21 ) 



,v 

Y j a k l<{‘~l<) = 0 1 £ / < N 



k=0 



where N is the order of the predictor being used and R is the autocorrelation function 
of the random sequence Y(n). The N equations generated by (21) can be written in 
matrix form for N = 2 as 



R-i Ro R - 1 

R 2 R i R 0 



1 


1 

1 




1 

o 






— 


o 




- a 2. 







( 22 ) 



The error variance c is given by 



c = = E[e n , }•(/;) - f'(«)] 

= Rle, v }'(//)] . 



(23) 



Substituting (18) into (23) yields 




}’(/»-*),}'(«) 



k=0 



(24) 



A' 



which in terms of the autocorrelation becomes c = 
this can be combined with (22) to yield 



For example, with N = 2, 



1 

o 


R - 1 


*- 2 " 


1 

c? 
1 




C 


R> 


R o 


R- 1 


«1 


= 


0 


A 




1 

0 


_ a 2_ 




. 0 . 



(25) 



This autocorrelation matrix is seen to be (N + 1) by (N + 1), Toeplitz, and since 
= R k , it is symmetric. The various solutions for (25) occupy much of the literature 
on linear prediction. Tor our purposes, the desired solution must provide a simple ex- 
tension to two-dimensional signals and be adaptable to a lattice structure. J.H. Justice 
(Ref. 2] has derived such a solution in terms of Szego polynomials. Although quite 
similar to the Levinson algorithm, it has the important advantage of not relying on the 
symmetry of the correlation matrix. The two-dimensional autocorrelation matrix will 
be found to be Toeplitz in blocks, and, although this provides some simplification, the 
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symmetry of (25) is lost. Fortunately, a solution in terms of orthogonal Szego 
polynomials is easily derived [Ref. 3]. 

3. Applications of the Autoregressive Model 

Linear prediction plays an integral role in the solution of a number of engi- 
neering problems. The resulting model can be used for high resolution spectral esti- 
mation [Ref. 4]. Speech signals can be reduced from a large volume of data to a 
comparatively small data space by linear predictive coding [Ref. 5]. A model of the 
generating function for the random sequence can be derived. Lattice filters can be im- 
plemented which can be used to form the prediction-error filter in the forward direction 
or synthesize the random process in the inverse direction. 

The many applications of linear prediction in time series analysis has generated 
a great deal of interest in extending the results to multi-dimensional random processes. 
Multi-dimensional signals occur in situations where the samples are spatially dependent 
such as radar, sonar, and image processing. The signals considered in this thesis are wide 
sense stationary and scalar-valued sequences of two variables. The two-dimensional 
autoregressive model that is derived will be used in some of the same applications clas- 
sically associated with the 1-D model including spectral estimation and sequence gener- 
ator estimation. An alternative implementation of the model in a lattice structure will 
also be presented. 

C. NOTATIONAL CONVENTIONS 

1. Matrix Notation 

This thesis will deal with a finite size rectangular sample space. The individual 
samples will be denoted by Y(n { —k, n 2 — /), where k and / are measures of distance from 
the point at which the value of the sequence is being estimated. The sample space will 
contain XI x N2-1 samples. Figure 2 shows an example of the sample space with 
\1 = N2 = 4. In this case we would be predicting the value of Y(3,3) based on the 15 
samples available. Notice that, in this case, k and / have a maximum value of 3. Spa- 
tial samples are indexed by nl and n2 to denote rows and columns, or the vertical and 
horizontal directions, respectively. 

2. Causal Samples 

In one-dimensional signals the clear meaning of past and future leads to a 
straightforward definition of causality. In two-dimensional signals the ordering of the 
points is of central importance; therefore a reference for past and future needs to be es- 
tablished. In this thesis causal samples will be taken to be in the quarter plane shown 
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Figure 2. Sample space for N 1 = N2 = 4. 



in Figure 3. This is equivalent to describing the desired filter as a third-quadrant filter. 
Due to the symmetry of the 2-D autocorrelation the third and first quadrants produce 
the same results. However, causal will be used only in reference to the samples shown 
in Figure 3. 

3. Polynomial Notation 

Two-dimensional polynomials will be denoted by P„, where k is the degree of the 
polynomial in the nl direction and / is the degree in the n2 direction. When the need to 
refer to an individual clement in a polynomial arises, the element will be identified as 
p kl (ij), where i and j give the power of the bivariate variable. The correspondence be- 
tween the z transform of the system model and the coefficients in the orthogonal 
polynomials leads to the use of z, and z 2 as the independent variables in this study. 
Thus the polynomial P 22 will be a 3x3 bivariate polynomial consisting of the elements 
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A 




data 



Y(n1,n2) 



Figure 3. Causal Data Space 



1*22 






7^22(0, 0 ) 



^ 22 ( 0 , l)r 2 



z - J /> 22 (2,0)z 1 2 



7’22(^'0 2 ! 2 2 

PllVWA 2 ! 



Pn(.Q’2)z 2 
P 22(*>2)z,r 2 
7-22(2, 2 ) 2 | 2 z 2 2 



(26) 



The development of P„ is in terms of positive powers of z. The desired filter 
must be causal and therefore in negative powers of z. This is easily accomplished by 
multiplying P*, by zr*Z 7 '. The resulting causal filter will be denoted A k , . The leading 
coefficient of P kl is p k ,{k,l ) which will become a^CfO) in the causal filter. 
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II. TWO-DIMENSIONAL LINEAR PREDICTION 



A. PROBLEM DEFINITION 

The two-dimensional linear prediction problem is conceptually identical to the one- 
dimensional problem. That is, we attempt to predict the value of Y(nl,n2) based on a 
linear combination of the previous NlxN2 rectangular block of data. The coefficients 
which result in the optimum prediction will comprise a linear, shift invariant, causal, 
prediction error filter. The problem is to find the coefficients that minimize in some 
sense the error between the actual value of Y(nl,n2) and the predicted value, F(nl,n2). 
This will require the solution of a block-Toeplitz matrix which contains NlxNl blocks 
and each block will contain N2xX2 correlation coefficients. We will require the solution 
to provide the means for a direct realization of the filter or implementation in a lattice 
structure. The mean square error criterion is again utilized to obtain the optimum pre- 
diction filter. If a solution to the desired filter can be achieved recursively, then the re- 
lationship between successively higher order filters can be used to implement the lattice 
filter. 

We begin by seeking a solution to the linear prediction problem, and will assume 
that Y(nl.n2) is a two-dimensional stationary random sequence. Given the causal block 

A 

of XlxX2 data samples we need to generate T(nl.n2) from a linear combination of the 
data that minimizes the prediction error. The prediction, f'(nl,n2), can be written as 

T(«l.«2) = - a 10 Y(nl - 1,«2) -c 01 («l,«2 - 1) - a n T(«l - 1,«2 — 1) - ••• 

- — — *'(«! - ' VI + E«2 “ N2 + I) 



where the a s ,'s are the filter coefficients associated with data samples that are delayed by 
.s in the nl direction and t in the n2 direction. The coefficients can be assembled in a 
matrix to yield the desired two-dimensional linear prediction filter. This filter will be an 
all-zero model with the finite sample space Y as the input and the prediction, Y, as the 
output as described below. 
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Y 2-D filter coefficients (28) 



— 1,A'2— 1 


a A' 1 — 1,1 


*AM 


a l,A r 2-l 


*11 


*10 


*0,A’2— 1 


*01 


0 




Figure 4. Generation of the prediction error 



The prediction process can be depicted as in Figure 4. The polynomial A(r, , z 2 ) is 
simply the z transform of (27) given by 

/J(2,,2 2 ) = - (fliozi -1 + *oi 2 2~' + a \i z \~' z 2 7 + -•) • (29) 
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For any choice of the coefficients some error will be present. This prediction error is as 
seen in Figure 4 and can be written as 

(VI —1 <V2— 1 

e = Y(n 1 m2) + EE a sl Y(nl - s,n2 - t) . (. s,t ) * (0,0) (30) 

j=o r=o 

We will define the filter coefficient a oa = 1 and condense (30) to the form 

A'l— 1 V2-1 

« = E Y J a ^ ni ~ s ’ n2 ~^ ■ < 31 ) 

j=0 f=0 

The task now is to find the filter coefficients that minimize the error. 

B. TWO-DIMENSIONAL ORTHOGONALITY 

1. Description 

The error, Y(nl,n2)-F(«l,«2), can be visualized as in Figure 5. Any linear 

A 

combination of Y(nl-l,n2) and Y(nl,n2-1) will yield an estimate of Y which lies in the 

A 

plane of Y(nl-l.n2) and Y(nl.n2-1). If Y(nl,n2) lies in this plane, then Y will be exact 
and there will be no error. For any other Y(nl,n2) the error will be minimum if it is 
orthogonal to the sample space spanned by Y(nl-l,n2) and Y(nl,n2-1). This is equiv- 

A 

alent to defining Y as that portion of Y(nl,n2) that is parallel to the sample space. Thus 
Y- Y leaves only the component of Y that is perpendicular to the sample space as the 
error. Although it can not be drawn for higher order predictors, the minimum error is 
always such that it is orthogonal to the space spanned by all the samples of Y used in 
the estimate. 

2. Derivation 

The orthogonality of the error to the sample space can be easily derived. First 
we define the inner product, <X,Y> , of two sequences, X and Y, by 

< x , y > =EE a »' 1 *'- (32) 

k 1 

Now the inner product of the squared error can be expressed as 

<e 2 > = <(F— }f> (33) 

and we seek a minimum squared error with respect to the filter coefficients a. 
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Figure 5. A visualization of the orthogonality principle for tno samples. 



~4r < e 2 > = 0 r (34) 

da 



or 



< 2e — rr e > =0 (35) 

da 

where —r are the partial derivatives with respect to the elements in the coefficient ma- 
da 

trix. Substituting (30) for e 




/VI -i 



-4- [>'(»!. w 

da 




A r 2 — 1 

Yj a st Y ( ni ~ 5 ’" 2 ~ 03 



t-Q 



(36) 



leads to 



13 



~^re = [ Y{n 1 - 1 ,n2),Y{n\,n2 - 1), •••• , Y{n\ - AT - 1,«2 - N2 - 1)] . (37) 



da 



Substituting (37) into (35) yields 



< e,Y{n\, — s,n2 — t)> =0 , 0 < s < AT 



0 < / < N2 
s,t=£ 0 . 



(38) 



This result, (38), states that the inner product of the error and each of the past 
samples of Y(nl,n2) used in the estimate is zero which requires that the error be 
orthogonal to the sample space. Not only is the error orthogonal to the sample space, 
but errors created by successively higher-order predictors are orthogonal to each other. 
Also, the sequence e(n t , n 2 ) will be uncorrelated to the degree that the order of the pre- 
dictor allows. That is, if we had access to an infinite causal data space, then we would 
produce a purely white error sequence. Since we will always be limited to some finite 
size filter, we will produce an error sequence that is partially correlated. 

3. Two-Dimensional Normal Equations 

The sequences considered in this thesis are taken to be wide-sense stationary. 
If they are also taken to be zero-mean signals, the inner product of two sequences is 
equivalent to the expected value of the product of the sequences. The expected value 
of Y(nl,n2) and Y(nl-s,n2-t) reduces to a two-dimensional autocorrelation that is solely 
a function of the delays s and t. Thus, 

< Y{n\,n2),Y{n\ - s,n2 - t)> = £[T(«1,«2), T(«l - s,n2 - t)] = R sl (39) 

where R will be used to denote the autocorrelation with lags s in the nl direction and t 
in the n2 direction. Equation (38) can now be used to form a block-Toeplitz matrix that 
is analogous to the Toeplitz matrix generated by the normal equations in one- 
dimensional linear prediction. Thus, we can write 



E[_e, }'(/?! — s,n2 — t)] = 0 (s.t) ^ (0,0) 



(40) 



substituting for e 



M-l :V2— i 



E F(wl,/?2)+ ^ ^ a kl Y{n\ — k, n2 — /),}’(«!— s, n.2 — t) 



= 0 0<s<Nl 



k = 0 1=0 



0 < t < A’2 (41) 

M # (0.0) 
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that can be written as 



A'l-l A'2-l 0<S<NI 

Y = 0 . 0 < r < A'2 (42) 

0 /=0 s,t * 0 

For example, with a N1 = N2=2, and with the First filter coefficient defined as 
a 00 = 1 the normal equations become 

E I Y{n 1 ,n2) + u 01 Y(n\,n2 - 1) + a 10 F(/jl - 1,«2) + 

+ a u Y(n\ - \,nl- 1) , F(«l -s,n2-t)\ = 0 

where (42) or (43) are the two dimensional normal equations. They are a direct result 
of the orthogonality relationship described by (40). These equations can be written out 
to yield NlxN’2-1 linearly independent equations. For example, with N1 = N2=2, we 
have 



s = 0 ,r = 1 


*01 + * 01*00 + * 10 *-,., + *11 *- 1.0 = 0 


(44) 


s = 1 ,/ = 0 


*10 + * 01 * 1, -1 + * 10*00 + * 11 * 0, -1 = 0 


(45) 


5 = l./= 1 


*11 *F * 01*10 "b * 10*01 *F * 11*00 = 0 


(46) 


Equations (44)-(46) can be augmented by the expression for the squared 
is given by 


error 




£ = £ 0 2 ] = E[e, Y- F] 


(47) 



or 



A'l-l A'2— 1 



= E 



y, a ki y ( n 1 - - /). Y(n 1 ,n2) 



(48) 



k = 0 1=0 



and in terms of the autocorellation 



A'l-l A'2— 1 

= E L a " Ri <■»> 

k = 0 1=0 



Assembling the above equations into a matrix yields 
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( 50 ) 



I 

T 

7 

CD 

7 

7_ 

cT 

o 

o 

1 




a 00 




£ 


d f 3 

c? 3 

0 

1 
i 

o 




^01 




0 


7 

o 

o 

7 

o 




a 10 




0 


1 

o 

o 

o 

o 

l 




_« 11 _ 




_0_ 



The resulting matrix is (NlN2)x(NlN2) square and will consist of blocks which 
are Toeplitz. There will be NlxNl blocks and each block will contain N2xN2 corre- 
lation coefficients. The matrix can be seen to be Toeplitz with respect to the blocks and 
can be written as 



®0 


of" 


% 




£ 


_ 4>1 


O 0 


_« 1 _ 




_ 0 _ 



(51) 



where the O/s are blocks in the autocorrelation matrix and H indicates hermitian 
transposition. Each O* is an N2xN2 Toeplitz matrix. For the example from (50), we 
have 



*00 


Vi 


*o> 


*00 


*10 


*1.-1 


*11 


*10 


(^00, 


«o.) 7 ' 


(«10, 


a \\) T 


: ( c . 


0 ) r 



(52) 

(53) 

(54) 

(55) 

(56) 



If an element in O* is R s , , the corresponding element in O w is R'_ Stl . Since we 
are dealing only with real data the conjugate notation can be dropped from R. It is also 
worthwhile to note that R s , = R _ s _, . 

The solution of this matrix for the coefficients a u will yield the two-dimensional 
whitening filter. There has been a great deal of research involving block Toeplitz ma- 
trices and there are several algorithms that could be used to obtain the coefficients. In 
order to implement the filter in a lattice structure, a recursive solution beginning with a 
zeroth- order filter and proceeding to the desired order is necessary. In one-dimensional 
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prediction error filtering, the Levinson algorithm successively represents the higher order 
filters in terms of reflection coefficients. It will be seen that an analogous solution to the 
two-dimensional prediction error filter can be developed in terms of Szego polynomials. 

C. SZEGO POLYNOMIALS 

The solution to the system of equations Oa = £ does not require that the inverse of 
the matrix <f> be found. All that is required is an orthogonal set of polynomials that span 
the space defined by O. This result follows the orthogonality relationship of the pre- 
diction error and the sample space. For example, if we had a second order and a first 



order prediction error filter the corresponding prediction errors would be 

c, = F(l)-«,(imO) (57) 

e 2 = F(2) — a 2 (l) E(l) — a 2 (2) F(0) (58) 

£[e 2 .c,] = E[e 2 . F(l) — a,(l)F(0)] (59) 

But (44) requires that e 2 is orthogonal to both Y(l) and Y(0), and therefore 

Ele 2 c,] = 0 . (60) 



It is evident from the above observations that the choice of A*., must be such that 
it is orthogonal to A*. If we begin by defining A 0 as any convenient value and proceed 
to find the successively higher-order orthogonal polynomials until the desired predictor 
order is reached, the final polynomial will hold the coefficients which will provide the 
optimum prediction-error filter. To accomplish this we will use Szego polynomials. 
Szego polynomials are a special class of polynomials which have zero inner products 
over an interval on the unit circle [Ref. 6]. This is precisely the requirement of the 
orthogonal polynomials relating successively higher order-prediction errors. J.H. Justice 
[Ref. 3] has presented a recursive method of generating Szego polynomials of the form 

p k = P(kf k + Pk-\ zk ~ X + -• + Po (61) 



where 



z = (62) 

The matrix O is an inner product matrix which may be utilized to orthonormalize a 
sequence of polynomials on the unit circle [Ref. 6]. The resulting orthonormal sequence 
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will be composed of Szego polynomials. If the filter sequence A* is determined in terms 
of Szego polynomials that are orthogonal over the space associated with the 
autocorrelation matrix , the prediction error filter will result [Ref. 2]. Further, the Szego 
polynomials can be determined in a recursive manner which is necessary for the repre- 
sentation of the filter in a lattice structure. In a one-dimensional filter this will result in 
the Levinson algorithm. J.H. Justice [Ref. 2] has developed the bivariate extension of 
Szego polynomials that can be utilized to ortho normalize the block Toeplitz matrix 
generated by the two-dimensional autocorrelation of the sequence Y(nl,n2). 

1. Levinson Algorithm in Terms of Szego Polynomials 

Any positive definite matrix can be defined as an inner product space and can 
thus be used to define orthogonal polynomials spanning that space. The autocorrelation 
matrix produces such a space. If the autocorrelation matrix is generated from a se- 
quence, 7= F 0 , 7,, • ,Y S , then the sequence of polynomials, 1 ,z,z 2 , ,z* can be 
orthonormalized with respect to the autocorrelation matrix. The resulting polynomials 
are P 0 (z), P^z), , Pp(z). The subscript on P denotes the degree of the polynomial. The 

requirement on P* is that it is orthogonal to all polynomials of lesser degree. Thus, we 
have 



A one-dimensional sequence of orthogonal polynomials will consist of the terms 



< P k , P,> = S(k-[) 



(63) 



Po~ Poo 

P\ =P\\Z + PlO 



(64) 




and in general, 



* 




( 65 ) 
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Justice's method for determining these polynomials recursively results in an al- 
gorithm equivalent to the Levinson algorithm. This algorithm for the 1-D case is briefly 
summarized as follows. We begin by setting 

P 0 = tf 0 "T . (66) 

Then the subsequent polynomials can be determined by the recursive relationship 

P„+,M - (1 - \>jV m )-^i?Pn^)->-nP m P'M (67) 



where 



P r n (z) = z n P n (z~ x ) 



(68) 



and 



'Ll ^ ' J ^k+]Pnk • (69) 

k= 0 

The equation (67) can be expressed in the form of a causal prediction-error filter 
that yields the familiar form of the Levinson algorithm. The nth-order polynomial 
produced by (67) will be of the form: 

P n {z) = p nn 2 n + Pnn-\Z n ~' + -• +PnO • ( 7 0) 

The causal prediction-error filter is of the form: 

^ n( z ) = a nO + a n\ : ^ + a nn z " ( 7 0 

where we let a np = p njx _ p and A„ = z~ n P n . Then multiplying (67) by z _(n+l> yields 

z~ {n+X) P n+x {z) = z~ n P n {z ) - z~ {n+l) ;. n P r n (z) (72) 

where the normalizing factor, (1 — U„l , has been absorbed into the definition 
of and the leading coefficients are defined as = a n0 = 1. Now, we define the 
backwards polynomials as follows 

P r n {z) = z n P n [z~ X ) (73) 
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(74) 



A n( z ) = z n A n (z ’) 



= z- fl (z n P„(z-')) 



p r n (z) = z n A r n {z) 

where (72) can now be written as 

A n +l( z ) = A n( z ) - 2 ~ l J.„ A '(z) 

A n+\( z ) = z ~ {n+ 1 ) { A n( z ~ l ) - zX H A r H (z~ 1 )) 
= z~'A r n (z) - z~ n ?. n (z n A n (z)) 

A n ( z ) = z ~' A n( z ) ~ K A A Z ) 



(75) 



(76) 

(77) 

(78) 



7 ,^, + i — /) 



4 = 



/=0 



/=0 



R 



(0 



(79) 



Equations (76), (78) and (79) represent the usual Levinson recursion of one- 
dimensional linear prediction problems. The coefficients found using the Levinson al- 
gorithm will be developed in negative powers of z. The coefficients found using (67) are 
in positive z and differ by a normalizing factor. Both sets of polynomials, A and P , will 
be orthogonal sets of polynomials on the unit circle spanning the same inner product 
space. Since the definition of the is simpler using (67), we will proceed to the two 
dimensional extension using positive powers of z and make the necessary conversion to 
a causal filter after the coefficients have been found. 

2. 2-D Extension 

Now we want to extend these ideas to two-dimensional filters. The necessary 
extension to bivariate orthogonal polynomials has been developed by Justice (Ref. 2). 
The derived polynomials will be orthogonal in two variables with respect to the inner 
product space of O as described by 



ElP kh P s J = 3(k-s)S(l-t) . 



(SO) 
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The set of polynomials to be orthonormalized are shown in Figure 6. The procedure 
starts with the first row which is a function of z 2 alone. This row, therefore, can be 
orthonormalized by the same method presented for one-dimensional systems. 







1 

z 2 
1 1 

zz 
1 2 
2 1 

ZZ 

1 2 



2 

z 2 
1 2 

^ z 2 

2 2 

ZZ 

1 2 



• • • • 



• • • • 



• III 















Figure 6. Variables used to form the bivariate orthogonal polynomials 



This is equivalent to taking the block <f> 0 from the two-dimensional autocorrelation ma- 
trix and finding the filter coefficients a 01 , <r 02 , - , c^ m . Thus, 



% 



*00 *0.-1 *0,-m 

*01 *00 •••' *0,-.(m-l) 



(81) 



*m 0 



*00 



J 



where 



a o ~ (°oo> «oi- •• > a Q m) 



(82) 



It is easily shown that 7^ 0m = R 0 _ m , so that this is a symmetric matrix and a re- 
cursive solution has been presented. Once the first row of orthogonal polynomials has 
been determined, the polynomials of the next row can be expressed as a linear combi- 
nation of the lower-order polynomials. For example, with N2 = 3 the first row will have 
polynomials P 00 , P 0I , and P 02 . The first polynomial in the second row to be found is 
P, o- This can be accomplished as follows: 
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First note that if P 00 is multiplied by z, its order is changed and it is no longer a 
member of the orthogonal set. We may then write 

z lA )0 = ^-oo^io + Joi^oi + ^ 02^02 (S 3 ) 

and after multiplying on both sides by P l0 and taking the expected value we find 

Woo, P loll = ;-00 W, P 10 ] + A 01 £[Poi^ 10 ] + ^-02^1^02, ^ 10 ] • (84) 

Since P I0 is orthogonal to all polynomials of lower order in z„ this reduces to 

EL Z \P oo, ^io] = ^oo (85) 

and similarly 

Woo. P 0i] = ^oi (86) 

Woo, A> 2 ] = ^02 • (87) 

Thus, in general we find that 



/-I k A7-1 

z i Pki = hlPi k+\,l + E 

r=0 s = 0 r=0 

Now we define the reverse polynomial 

p;p„z 2 ) = zf^zrWVo^'s* 

0 <j < m (89) 

m = N2 — 1 

Note that P r tJ will also be an orthogonal set of polynomials and will span the same space 
as P iy Now (8S) may be written in terms of the reverse polynomials 

/-I k m 

z \Pkl = ^kl^k+XJ + ^/kiP fe+l,f + A (90) 

/=0 s =0 r = 0 

where 

r„ = E[z,p t „ p;j . (9i) 

But for 0 < s < k — 1 



ktP k+ \ ,t + 



E E a 



( 88 ) 



~n 



£[z,p«.' p ;,]= o 



(92) 



so that (90) may be reduced to 



z \Pkl~ '‘■kl^k+\,l + 






r= 0 



m 

t = 0 



(93) 



where 



— E[ z \Pkb Pk+i,tl ( 94 ) 

and 

2'*, • (95) 

The polynomials, P MJ , may now be found from the lower order polynomials 



/-i 

VW = z t^w — 7.VW — 

2=0 2=0 



(96) 



ft+i,/ 



^■kl^k+lj 

II^A/^+1,/1 



(97) 



This completes the development of the bivariate Szego polynomials as pre- 
sented by Justice [Ref. 2]. The polynomial P*, will hold the coefficients for the k xl 
prediction error filter. The desired causal prediction error filter, A*„ is obtained by 
multiplying P w by zy^zy' and dividing by p k ,{k y [) to force the leading coefficient 
a*,{0,0) = 1 as was assumed. The signal generator model is found from (30) with efortj) 
taken as the input. 



*(m,z 2 ) 



1 

^( 2 i> : i) 



l_ 

+ « 10-1 1 + a 01 2 2 1 + a 



11 



-1 -1 
■1 2 2 



+ 



This again assumes that A(z,, z 2 ) is a causally invertible sequence. 



(98) 
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D. LATTICE IMPLEMENTATION 



1. Analysis Lattice 

It is now a simple matter to utilize the orthogonal polynomials for lattice im- 
plementation. The first row is a function of z 2> only and the lattice structure is identical 
to the one-dimensional case. All subsequent polynomials have a connection with each 
of the polynomials in the previous row and each of the polynomials in the same row but 
of lower order in z 2 . 

Using the familiar form of the Levinson algorithm the first row of the lattice 
filter may be drawn as shown in Figure 7, given by 

^O./t+l = ^0 k ~ z 2 ^-k^Ok (99) 

and 

^0,/t+l ~ z 2 '^o k ~ S-k^ok • (100) 



The orthogonal polynomials that have been developed will hold the filter coefficients in 
descending powers of z, and z 2 . That is, the lowest order term in P k , is the term associated 
with the greatest delay in the filter difference equation. We could use (94), (95) and (96) 
to implement a first-quadrant filter that would be statistically equivalent to the third- 
quadrant causal filter we desire. If all the data are stored ahead of time, which is gen- 
erally the case for 2-D sequences, this would not cause any computational difficulties. 
We will present the causal implementation of the lattice structure. The first quadrant 
development is similar. 

In order to correlate the powers of z, and z 2 with the delay and thus have a 
causal filter, we will define A k , — z\ k z 2 'P kl . This will result in a realizable design and the 
form of the lattice filter. The first row of the filter is given by (99) and (100). We now 
assume that row k has been constructed and that we proceed to row k+ 1. Multiplying 
(96) by zf ( *‘ 1) z 2 *' yields 






(*+,)__/ 

kh 1 z 2 1 k+\,l ~ 



l - 1 



-1 L 2 1 kl 



-z 



; 7 ~ l p 

A kr\ ^2 1 



1=0 



- ^VZ-'-'iAaoi) 

r=0 



that can be written in the form 
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Figure 7. First row of the analysis lattice 
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^■kl^k+\,l ~ A hi 



l - 1 

-l- 



i=0 



a _ ,-»W r r 

/ -fcf /l /!+l.r 2 1 / / • kt A kt 

i=0 



(102) 



Now with A' uu defined as zf (A+1) Z 2 m A i+1 /zf', z 2 ‘), we have 



^■k/^k+UI ~ z i '^k 



k! 



I - 1 

2 

r=0 



4'-%, 



/ Wl,r 



rn 

-E 



2V/1 



ki^ki 



(103) 



r=0 



Equations (102) and (103) provide the means for implementing the two-dimensional 
lattice structure. Before proceeding to a schematic representation of the structure, it is 
worthwhile to explore the relationship of the 2*/s and the filter coefficients. This will 
provide some insight into the operation of the lattice filter and make it apparent that the 
2's are sufficient to represent the filter. Assume that we have already found the lower 
order filters, A 00 , A ou A 02 , and A I0 and now seek A u . The coefficients that are available 
are shown in Figure 8. By (102) we find that 




Figure 8. The lower order polynomials for A n 
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(104) 



''•01^11 — ^01 — z 2 1/J -00^10 — Z 1 ^-'oo^oo — Z 1 ^-'oi^Ol — Z 1 02-^02 



The filter A n (z v z 2 ) will then be of the form 



^ll( z l' z 2) — «00 + «01 z l 1 + «10 z 2 ■f a ll z l z 2 



(105) 



and (104) can now be used to equate the coefficients in A n with the lower order 
polynomials. This is shown in Figure 9. The equations for the coefficients can be writ- 
ten directly from Figure 9 as 



flu (00) = o 0 i(°°) 


(106) 


On(01) = 00,(01) - ;. 00 o 10 (00) 


(107) 


a ll(10) = — /i - , 02 a 02(02) 


(108) 


— E>-oo a io(10) + z' 0] a 0 ,(01) + )’ 02<^ 2 (01)H 


(109) 


^-’ooflooCOO) + / - , oi a oi(00) + ^-’oz^COO)] • 


(110) 



It can be seen that there are enough equations to determine all the filter coefficients of 
A„. Conversely, if the filters are known, then the Z.'s can be completely specified. This 
relationship allows the filter to be realized in either the direct form using the filter coef- 
ficients or in the lattice structure using the /'s. If the same procedure is carried out 
starting with (102). the coefficients in A r n are found to correspond to those already found 
for A n . This is of course a requirement if the algorithm is correct. 

The second row of the lattice structure described by (102) and (103) is shown in 
Figure 10 for a 2x3 filter. The first row of the filter is shown in Figure 7. All subsequent 
rows of the filter will have the same configuration. The structure can be seen to be quite 
similar to a one-dimensional lattice. There are. however, two important differences. 
First, there is a connection to all the lower-order predictors in the row. Also the inverse 
delay in the reverse error path appears to be a non-causal operator. These inverse delays 
result from the definition of A\,= zf*zj m /fi / (zf I , zj 1 ) and force a specific ordering of points 
being inputed to the lattice. For example, if we were computing the reverse error, 
e r u (n\,n2), for a filter with three horizontal filter coefficients, m = 2, then the z transforms 
of of these errors are given by 

£[,(zl,z2) = Y(z l ,z 2 )A[ l (z l z 2 ) (111) 
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( 112 ) 




£[ 0 (zl,z2) = F(z,, z 2 )(a 00 z, 'z 2 2 + a 10 z 2 2 ) . 



(113) 



By referring to Figure 9 we see that the coefficients a„(0,l) and a 10 (O,O) must be 
in the same power of z,z 2 and likewise a,,(l,l) and ai 0 (l,O) . Therefore, /![,(z,, z 2 ) must be 
multiplied by z 2 This would cause no problems if the direct form implementation was 
being used since 



There are no non-causal terms. In the lattice structure instead of linking actual filters 
we link error outputs. The calculation of e' u {n\,n2) requires e' m {n\,n2 + 1) . This requires 
the next horizontal value of e l0 before the present value of e n can be computed. The end 
result of this is that the previous row errors must be completely calculated and stored 
before the present row errors can be computed. This forces the input to proceed across 
rows. 

The filter shown in Figure 10 represents extending the filter order by rows. If 
at some point we wish to extend the filter by columns the procedure is identical and the 
next column would have the same connections as the row extension with the subscripts 
transposed. 

2. Synthesis Lattice 

The all-pole model provides the means to synthesize the original signal from the 
error sequence. The synthesis filter can also be implemented in a lattice structure. 
Multiplying the recursion relationships (102) and (103) by Y(z) we arrive at the recursion 
for the error sequences. 



z 2 £[ 0 (zl,z2) = y^z^tfooz, 'z 2 1 + a w z 2 ') . 



(114) 





( 116 ) 
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Figure 9. Coefficient relationship between orthogonal polynomials. 



e k+\,l 



-1 

h 



e kl 



/- 1 
/=0 




(117) 



Since Y(nl,n2) corresponds to e 00 , we must write (116) in decreasing order re- 
sulting in 



e k,l = e k+\,l + 




(118) 



The first row of this filter is again a function of one variable and is given by 

e kt ~ e k,i + 1 + ' > -i e ki z 2 1 



(119) 



or 

e k,l + 1 = e kl - }-l c kl • 



( 120 ) 
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Figure 10. One cell of a 2-D analysis lattice 
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This is realized as shown in Figure 11, and the subsequent rows are given by 
(118). The signal flow graph is similar to Figure 10 with the signal flow inverted and the 
sign changes as indicated by (118). The node at e 00 is shown in Figure 12. Computer 
simulation of these filters and the results obtained are presented in Chapter III. 
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Figure 11. First row of the synthesis lattice 
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^°° e(n 1r n 2 ) 
oo 



e(n.,n 2 ) 

01 



e(n 

0m 



Figure 12. node of the synthesis lattice 
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III. COMPUTER IMPLEMENTATION AND RESULTS 



A. IMPLEMENTATION 

1. Autocorrelation Function 

The autocorrelation function referred to throughout this thesis is actually an 
estimate based on a finite sample space. There are several ways to compute the esti- 
mated autocorrelation coefficients. For a large, complex-valued sample space the most 
efficient method is the use of Fourier transforms. For the small real-valued data space 
used in the examples of this study, the autocorrelation is estimated directly from the 
available data. The autocorrelation coefficients are given by 



2. Reflection Coefficients 

The basic algorithms for calculating the coefficients were presented by Justice 
[Ref. 2]. With some modification these algorithms are presented in this and the follow- 
ing section. 

The coefficients that link succesive order filters are derived in Chapter II. They 
will be referred to as reflection coefficients even though this term is descriptive of only 
the first row. The first row reflection coefficients were found to be 



This result is simply the one dimensional Levinson algorithm appended with a second 
subscript. 

For subsequent rows there are two reflection coefficients, one for the forward 
polynomials and one for the reverse. They are obtained as 




( 121 ) 



k 




( 122 ) 



>-ki= < z\ p kb p k+\,t > 



(123) 



and 




(124) 
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These are inner products over the autocorrelation space and are computed as 



1 k + 1 r 

h, - k+\-rJ-sPk+\,t( r ’ (125) 

j=t r =0 5=0 

and 

m k I 

*'kt = Tsk^ Yj R -( r +')j-s pk ^ r ’ s ^ ( 1 26) 

j= 0 r= 0 J=0 

where any negative value in a subscript results in a zero. 

3. Filter Coefficients 

The First row of prediction error filters are polynomials of the form P 0k . They 
can be calculated by the relationships 

f 0 ,o(0,O) = (/J 00 )-T (127) 

/>(),/+ i(°J) = ( ] - I >•/! - 1) - '‘■iPo,i(°APoP’ 1 -j)) ( 12S ) 



0 < j < l + 1 
0 < / < m — 1 



( 129 ) 



where m is the maximum filter order in the n2 direction and l is the filter order presently 
being computed. 

All subsequent rows will increase the filter order in the nl direction and can be 
recursively computed by 



/-i 



m 



Pk+u(‘xi) = Pki‘ - it I) - /_J-k<Pk+ i jf (v) - A r kiPkl k - '> -j) 






i—m—j 



( 130 ) 



0 <i<k + 1 

0 <j < l 



( 131 ) 



This recursion can be carried on until the desired filter order is reached. The only 
polynomials required to increase the filter size are the preceding row and the lower order 
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filters in the same row. It may be desired to normalize the filter equation, and this is 
easily carried out by 




pk+U'J) 



(132) 



where 



<4+i,/ ~ II ^4+i, /II 



(133) 



or 




•*+i / 



(134) 



B. RESULTS 

1. System Model 

The algorithm which has been developed in this study is particularly useful in 
finding the autoregressive or all-pole model of the system which has generated these 
data. Several examples are presented in this section. Notice that the computer gener- 
ated results have the negative of the desired values. This results from our original defi- 
nition of the prediction of }' in (9) as 



Figures 13 and 14 show various size system models that were derived from the 
data generated by the difference equation 

Y{nl,n2)-.4Y(nl - M2) + .6F(/j1,«2- 1) - .6Y{n\ - M2- 1) + A'(«M2) . (136) 

Figure 13 contains the system models resulting from driving the difference equation with 
an unit impulse. Figure 14 is the result of driving the system with a pseudo white noise 
sequence. The discrepancies in the filter coefficients can be mainly attributed to insuf- 
ficient 'whiteness' in the noise data. 

2. Spectral Estimation 

If we assume the signal, Y(nl,n2), is generated by the model of Figure 15, the 
power spectrum is given by 



Y= - [a 01 Y{ti 1 ,w2 - 1) + tf 10 F(«l - M2) • • ] . 



(135) 
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FILTER ORDER— N1=2,N2 = 1 
1.0000 -.6000 

-.4000 0.6000 

0.0000 0.0000 



FILTER ORDER— N1=2,N2 = 2 



0 .0000 
0 .0000 
0 .0000 



FILTER ORDER— N1=3,N2 = 0 
1 .0000 
-.0917 
0 .0000 
0 .0000 



1 . uooo 
- . 4000 
0 .0000 



-.6000 
0.6000 
0 . 0000 



-.0001 

-.0000 

0.0001 



FILTER ORDER— N 1=2, N2 = 3 
1.0000 -.6000 0.0000 

“•4000 0.6000 0.0000 

0.0000 0.0000 0.0001 



FILTER 


ORDER — N1=3,N2=1 


1 .0000 


-.6000 


-.4000 


0.6000 


-.0000 


0 .0000 


0 . 0000 


-.0000 



FI 


LTER 


ORDER— N1=3,N2 = 2 




1 . 


0000 


-.6000 


.0001 


— , 


4000 


0.6000 


.0000 


— , 


0000 


0.0000 0 


.0001 


0 . 


0000 


-.0000 


.0000 



FILTER 


ORDER — N1=3,N2 = 3 




1 . 0000 


-.6000 0.0000 


0.0000 


- . 4000 


0.6000 0.0000 


0.0000 


0.0000 


0.0000 0.0001 


0.0000 


0.0000 


-.0000 0.0000 


0 . 0000 



Figure 13. System models for Eq.136 with an impulse input 
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FILTER ORDER — N1=2.N2=1 

1.0000 -.5407 

-.4136 0.5908 

0.0107 -.0274 



FILTER ORDER — N1=2,N2 = 2 

1.0000 -.5572 0.0111 

-.4125 0.6135 -.0404 

0.0099 -.0327 0.0130 



FILTER ORDER— N1 = 2,N2 = 3 

1.0000 -.5574 0.0070 

-.4349 0.6126 0.0007 

0.0284 -.0398 -.0302 



0.0307 
-.0287 
0 . 0699 



FILTER ORDER — N1=3,N2 = 0 

1 .0000 

-.2071 

0 . 0564 

-.0402 



FILTER 


ORDER-- N 1=3, N2 = l 


1 . 0000 


-.5402 


-.4132 


0.5903 


0 . 0074 


-.0216 


-.0043 


-.0133 



FILTER 


ORDER — H1 = 3,N2 = 2 


1 .0000 


-.5571 0.0139 


-.4126 


0.6129 -.0406 


0.0073 


-.0206 0.0082 


-.0052 


-.0188 0.0113 



FILTER 


ORDER— N1 = 3,N2 = 3 




1 .0000 


-.5570 0.0082 


0 .0302 


-.4365 


0.6125 0.0012 


-.0285 


0.0283 


-.0300 -.0469 


0 . 0781 


-.0130 


-.0126 0.0227 


-.0190 



Figure 14. System models for Eq.136 with a white noise sequence input 
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S(z v z 2 ) = (FlY(nl,n2)l) 2 



(137) 



or 



S(z u z 2 ) = (F\_B{n\,n2) * e{nl,n2)l) 2 



(138) 



or 



S(z,, z 2 ) = (2?(z t , z 2 )£(z,, z 2 )) 2 



(139) 



If the input is a unit variance white noise sequence, the power spectrum be- 
comes 



We can now use the all-pole system model to arrive at an estimate of the power spec- 
trum 



where A(z„z 2 ) is the all-zero prediction error filter previously derived. 

We have considered three examples of two-dimensional all-pole power spectrum 
estimation. In all examples the data set size used to compute the autocorrelation func- 
tion is 32 x 32. In the first example a sequence consisting of two cosines at normalized 
frequencies (0.5,0.25) and (0.75,0.75) is corrupted by two-dimensional unit variance 
white noise. Figure 16 shows the surface plot of the estimated power spectrum with a 
4x4 filter mask. The ideal response would be two impulse functions at the appropriate 
frequencies. The plot of Figure 16 shows several spurious peaks along a ridge which is 
illustrated by the contour plot in Figure 17. It can clearly be seen from the surface and 
contour plots in Figures 18 and 19 that the estimation of the power spectrum improves 
rapidly with an increased filter size. Figures 20 and 21 show the results using an 8 x 8 
filter. The surface plot in Figure 20 clearly shows two discrete peaks in the power 
spectrum. Notice that there is a slight bias in the estimates. 

In the second example three cosines at frequencies (0.5,0.25), (0.75,0.75), and 
(0.25,0.75) are used to generate the data. The surface and contour plots of spectral es- 
timates are presented in Figures 22 through 27. The filter masks are square with sizes 
4 x 4, 6 x 6, and 8x8. The ability of these filters to distiguish the spectral peaks again 



•S(z,,z 2 ) = B 2 (z u z 2 ) . 



(140) 




(141) 
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improves as the filter size increases from 4 x 4 to 8 x 8. While the 4x4 filter, as shown 
in Figure 22, provides very little useful information, the 8x8 filter in Figure 26 clearly 
indicates three distinct peaks. 

The third example investigates the capability of the algorithm to differentiate 
between two closely spaced sinusoids at frequencies (0.4375,0.4375) and (0.5625,0.5625). 
As observed in Figures 28 through 31 the algorithm has not been able to distinguish 
between the two peaks. However, with a filter size 8x8 (see Figures 32 and 33), the 
spectral peaks become distinct even though there are some weak spurious peaks present. 

The examples which are presented in Figures 16 to 33 show the usefulness of the 
algorithm in estimation of the power spectrum. The ability of the algorithm to identify 
discrete peaks in the power density spectrum increases rapidly with filter order. As can 
be seen in Figures 32 and 33, even closely spaced peaks are readily identified by an 
eighth-order filter. This filter requires the computation of 64 coefficients which is much 
simpler than using Fourier transforms. All plots are scaled from 0 to 1 which represents 
the interval 0,7r. 

3. Signal Synthesis 

The reflection coefficients for a single frequency data set were determined and 
used to implement a computer simulation of the synthesis lattice. The resulting lattice 
was then driven by a pseudo random white noise sequence. The spectrum of the output 
from the lattice is shown in Figure 34 for a 4 x 4 structure. Considering the small size 
of the filter and the poor quality of the input noise sequence, the results are acceptable. 
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4x4 FILTER 




Figure 16. Surface plot of spectral estimate: Y(nl,n2)= cos(nl nil + n2 rt/4) + 

co$(nl n 3'4 + n2 n 3,4) 
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SIXV z<t> 



o 



4x4 FILTER 




Figure 17. Contour plot of spectral estimate: Y(nl,n2)= cos(nl n',7 + n2 7 r/4) 

+ cos(nl 7x 3/4 + n2 n 3/4) 
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6x6 FILTER 




0 O 



Figure 18. Surface plot of spectral estimate: Y(nl,n2)= cos(nl n,'2 + n2 n,4) + 

cos(nl n 3, 4 + n2 n 3/4) 
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$2 AXIS 



6x6 FILTER 



o 




Figure 19. Contour plot of spectral estimate: Y(nl,n2)= cos(nl njl 4- n2 

+ cos(nl n 3/4 + n2 n 3/4) 
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8x8 FILTER 




0 O 



Figure 20. Surface plot of spectral estimate: Y(nl,n2)= cos(nl nil + n2 rc/4) + 

cos(nl n 3/4 + n2 n 3/4) 
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<t>2 AXIS 



8x8 FILTER 



o 




Figure 21. Contour plot of spectral estimate: Y(nl,n2) = cos(nl n\l + n2 ni 4) 

+ cos(nl n 3/4 + n2 n 3/4) 
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4x4 FILTER 




O 0 



Figure 22. Surface plot of spectral estimate: Y(nl,n2)= cos(nl nil + n2 ni 4) + 

cos(nl n 3, 4 + n2 n 3,'4) + cos(nl 7t/4, n2 n 3/4) 
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O . i ■ i ] 1 I i I . - I . I y I 

0 0.2 0.4 0.6 0.8 1.0 



*1 AXIS 



Figure 23. Contour plot of spectral estimate: Y(nl,n2) = cos(nl n'2 + n2 tt/4) 

+ cos(nl n 3'4 + n2 n 3.'4) + cos(nl tt/4, n2 n 3/4) 
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6x6 FILTER 




0 O 



Figure 24. Surface plot of spectral estimate: Y(nl,n2)= cos(nl nil + n2 7t/4) 4- 

cos(nl n 3/4 + n2 n 3/4) + cos(nl n! 4, n2 n 3/4) 
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$2 AXIS 



o 



6x6 FILTER 




Figure 25. Contour plot of spectral estimate: Y(nl,n2)= cos(nl n,'2 + n2 

+ cosfnl n 3 4 + n2 n 3 4) + cos(nI n’4, n2 n 3/4) 
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MAG 



8x8 FILTER 




0 0 



Figure 26. Surface plot of spectral estimate: Y(nl,n2) = cos(nl re/ 2 + n2 n,'4) + 

cos(nl re 3/4 + n2 re 3/4) + cos(nl re/4, n2 n 3/4) 
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$2 AXIS 



o 



8x8 FILTER 




Figure 27. Contour plot of spectral estimate: Y(nl,n2)= cosfnl n'2 + n2 n/4) 

+ cos(nl n 3,4 + n2 rc 3/4) + cos(nl n/4, n2 n 3/4) 
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4x4 FILTER 



16 



Figure 28. Surface plot of spectral estimate: Y(nl,n2)= cos(nl *7/16 -F n2 
7 t 7/ 16) + cos(nl*9/16 4- n2 *9/16) 
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$2 AXIS 



4x4 FILTER 



O 




Figure 29. Contour plot of spectral estimate: Y(nl,n2)= cos(nl 7x7/ 1 6 -I- n2 

7x7; 16) + cos(nl7i9, 16 4- n2 n 9/16) 
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6x6 FILTER 




Figure 30. Surface plot of spectral estimate: Y(nl,n2)= cos(nl 7:7/16 + n2 
7:7/16) + cos(nl7:9 : 16 4- n2 ti 9/16) 
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$2 AXIS 



6x6 FILTER 



o 




Figure 31. Contour plot of spectral estimate: Y(nl,n2) = cos(nl 7i7/16 + n2 

nl, 16) + cos(nlrt9'16 + n2 n 9/16) 
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8x8 FILTER 




Figure 32. Surface plot of spectral estimate: Y(nl,n2)=cos(nl *7/16 + n2 
7x7, 1 6) + cos(nl7t9/I6 + n2 7:9,16) 
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<t>2 AXIS 



8x8 FILTER 



2 1 1 1 1 1 1 1 1 r 
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Figure 33. Contour plot of spectral estimate: Y(nl,n2) = cos(nl tt7/16 

tt 7/16) + cos(nl7i9;16 + n2 n 9/ 16) 



+ n2 
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SYNTHESIZED SIGNAL 




Figure 34. Spectrum of signal generated from white noise: Y(nl,n2) = cos(nl n/2 

+ n2 tt'4) 
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IV. CONCLUSIONS 



The linear prediction approach to analysis of stationary one-dimensional signals has 
direct applications in the study of two-dimensional signals. The solution of the 2-D 
normal equations in terms of orthogonal polynomials provides an efficient and effective 
method for system modeling and power spectral density estimation. There also exist 
two-dimensional lattice structures that arc analogous to their one-dimensional counter- 
parts. These lattices can be used as an alternate implementation of prediction error fil- 
ters and signal syntesis filters. 

The simple nature of the recursive solution to the 2-D normal equations developed 
in this thesis provides an excellent method for obtaining solutions to problems that are 
not dependent on the stability of the derived filter. For example, system modeling 
achieves very accurate results assuming the system being modeled is stable. Signal syn- 
thesis on the other hand can not be considered reliable since there is nothing to guar- 
antee a stable filter. This can be attributed to the lack of autocorrelation matching in 
two-dimensional signals. An excellent discussion of this property can be found in [Ref. 
4J. Simply stated, there arc more independent autocorrelation coefficients than there are 
coefficients in the filters that we have found. We have ignored the second- and fourth- 
quadrant coefficients which arc independent from the first- and second-quadrant 
coefliccints we have been using. 

Considering the computational expense of other spectrum estimation techniques the 
algorithm presented in this paper is valuable if used only for that purpose. The problem 
of signal synthesis, while more difficult, was met with some success as evidenced by 
Figure 34. This limited success invites further study into autoregressive modeling of 2-D 
signals. 
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